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A.  Turan,  J.  Swithenbank,  S.A.  Billings  and  P.  Salib. 
Sheffield  University 


ABSTRACT 

Combustor  modelling  methods  involving  stirred  reactor  networks  require 
information  on  the  volumes,  inter-connections,  mixedness  and  flowrates  cf  the 
individual  reactors.  It  is  here  proposed  that  the  response  of  the  system  to  a 
pseudo-random  binary  sequence  (PRBS)  stimulus  can  be  used  advantageously  to 
fulfill  this  need.  A  survey  is  presented  of  the  theory  of  the  technique  and 
method  of  interpretation  of  the  results.  The  practicability  of  the  method 
has  been  investigated  experimentally  using  a  water  model  of  a  ga  turbine 
combustor;  salt  solution  tracer  was  introduced  into  the  primary  flow  in  a 
PRBS  and  the  conductivity  at  different  points  throughout  the  chamber  was 
used  to  determine  the  impulse  response  and  Z-transform,  The  results  confirmed 
the  anticipated  sensitivity  of  the  technique.  Finally,  the  extension  of  the 
technique  to  the  quantitative  identification  of  non-linear  elements  is 
presented. 


INTRODUCTION 


The  ultimate  objective  of  combustion  research  is  the  optimization  of 

practical  combustor  designs.  Typical  combustors  (e.g.  gas  turbine  combustors) 

involve  the  simultaneous  interacting  processes  of  three-dimensional  turbulent  flow, 

two-phase  evaporating  droplets,  mixing,  radiation  and  chemical  kinetics.  At  the 
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present  time  numerical  prediction  algorithms  are  becoming  available  ’  which  can 
model  all  these  processes  to  compute  the  hydrodynamic,  thermodynamic  and  chemical 
quantities  throughout  a  three-dimensional  field.  Complementary  stirred  reactor 
network  algorithms  permit  the  prediction  of  minor  constituents  (pollutants) , 
again  including  such  effects  as  droplet  evaporation  and  unmixednessz ’  ’  . 

Current  pressures  to  minimize  pollutant  formation,  improve  efficiency, 
tailor  the  pattern  factor  and  widen  the  turndown  ratio,  together  with  the 
developing  requirement  to  use  mpre  aromatic  fuels,  and  even  synthetic  or  coal 
derived  fuels,  suggest  that  the  initial  stages  of  combustor  design  would  be 
greatly  facilitated  by  reliable  mathematical  models  based  on  fundamental 
principles. 

In  order  to  represent  a  gas  turbine  combustor,  the  present  numerical 
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finite  difference  algorithms,  Fig.l.  typically  require  about  3000  grid  cells  , 
whereas  the  reactor  network  models  utilize  about  7  reactors.  Fig. 2. 

The  former  model  yields  details  of  the  spatial  distribution  but  would  require 
excessive  computer  time  if  complex  kinetic  equations  were  included.  The 
latter  model  forfeits  spatial  resolution  but  permits  complex  chemical  kinetic 
schemes  to  be  incorporated.  It  is  anticipated  that  the  models  can  be  eventually 
united  by  invoking  the  fact  that  optimum  local  grid  size  (for  both  models)  is 
largely  determined  by  the  turbulence  macro-scale.  Thus  spatial  inhomogencities 
in  the  fuel  feed  are  smoothed  in  a  single  stage  stirred  reactor  having  transverse 
dimensions  of  the  order  of  the  macro-scale.  Similarly  in  the  finite  difference 
formulation,  at  least  one  grid  point  is  required  in  a  shear  layer  whose 
dimensions  are  again  determined  by  the  turbulence  macro-scale. 


On  the  other  hand,  the  turbulence  micro-scale  determines  the  degree  of 

mixing  at  the  molecular  level  and  hence  the  maximum  possible  completeness  of 

combustion.  In  both  models,  the  micro-mixing  is  often  characterized  by  a 

turbulence  dissipation  parameter  such  as  the  eddy  break-up  paraneter^,  or 
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the  TgD  parameter  .  Whilst  the  t  parameter  will  easily  permit  the  prediction 
of  pollutant  formation  due  to  unmixed  and  hence  unbumed  material  by: 

\  •  >'<l  '  >'tsd>  »  “  -  l'tSD)  <» 

it  does  not  directly  allow  for  non-linear  effects  in  the  chemical  reactions. 

The  significance  of  these  effects  depends  on  the  amplitude  of  the  fluctuations 
and  the  extent  of  the  non-linearities.  In  the  case  of  N0X  formation,  it  is 
well  known  that  the  kinetic  equations  have  high  activation  energies  and  are 
therefore  very  non-linear.  Previous  efforts  to  account  for  this  effect^  have 
used  the  concentration  variance  g  to  compute  the  yield  of  N0x,  however  the 
problem  remains  to  predict  g  for  a  given  combustor  design. 

A  further  problem  arises  if  we  wish  to  characterize  a  single  reactor  in 
a  reactor  network.  Frequently,  there  is  strong  backmixing  or  re-cycle  in  the 
network  and  the  boundaries  of  a  single  reactor  are  often  drawn  round  the  whole 
volume  incorporating  the  backmixing.  It  is  here  proposed  that  the  longitudinal 
extent  of  a  single  reactor  is  determined  by  either;  the  distance  before  the 
stream  is  divided  (or  joined  by  a  second  stream),  or,  the  distance  by  which  the 
turbulence  decays  to  the  level  of  negligible  mixing.  The  transverse  extent  of 
each  reactor  is  determined  by  the  turbulence  macro-scale  rs  discussed  above. 
Backmixing  is  thus  defined  by  the  flow  rates  and  interconnections  between  the 
elementary  reactors,  and  parameters  such  as  blow-off  only  have  meaning  for  that 
part  of  a  network  incorporating  all  recycle  paths.  Each  elementary  reactor 
therefore  experiences  a  form  of  plug-flow  with  a  specified  degree  of  mi c To¬ 
rn  xing,  and  the  mixing  becomes  zero  in  the  conventional  plug  flow  part  of  the 
combustor,  thus  limiting  the  combustion  efficiency  as  indicated  in  Eq.V. . 

In  order  to  approach  this  ’ideal'  reactor  network,  experimentally,  powerful 
techniques  are  required  and  the  purpose  of  this  paper  i,s  to  present  some 
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appropriate  procedures  based  on  tracer  response.  The  method  is  also  > 

applicable  to  the  coalescence/dispersion  (C/D)  micromixing  mode 1 ^ . 

In  the  C/D  model,  micromixing  is  simulated  by  considering  the  reacting 

mass  to  be  made  up  of  a  large  number  of  equally  sized  'eddies'  behaving  as 

small  batch  reactors.  Eddies  undergo  collision  from  time  to  time  (at  a  mixing 

frequency  B)  and  equalize  concentrations  when  this  happens.  The  dimensionless 

parameter  I  *  8/a,  where  a  =*  l/(mean  residence  time  T  )  ,  is  analogous  to  T„„ 
m  s  SD 

as  shown  in  Fig. 3.  Computation  proceeds  by  randomly  selecting  a  pair  of  turbules 

within  the  reactor  which  average  their  properties  at  each  interaction  time 

interval  (t  /I  N),  where  N  is  the  total  number  of  turbules.  Whenever  a  feedtime 
s  m 

interval  (tg/N)  has  elapsed,  one  turbule  is  selected  at  random  to  be  removed 
from  the  reactor,  and  is  replaced  by  fresh  feed.  However,  as  indicated  by  Pratt, ^ 
provided  8  is  known,  the  C/D  model  can  also  be  used  to  obtain  both  the  mean 
and  standard  deviation  of  parameters  such  as  temperature  and  concentration  at 
the  combustor  exit,  together  with  their  various  correlations.  The  model  thus 
allows  for  non-linear  chemical  kinetic  rates.  In  the  case  of  the  TgD  model, 
the  conventional  treatment  can  be  regarded  as  reducing  the  reactor  volume  to 
allow  unmixed  material  to  pass  through  unreacted^ with  the  remainder  of  the 
reactor  acting  as  a  perfectly  3tirred  system.  This  is  reasonably  true  for  gas 
turbine  primary  zones  3ince  is  so  high.  If  the  variation  of  the  temperature 
and  concentration  about  the  mean  could  be  related  directly  to  T__,  then  by  analogy 

Ov 

with  a  8 /<!  model,  it  is  possible  to  divide  the  reactor  into  'cells'  operating  at 
appropriate  conditions  and  thu3  allow  for  the  non-linear  chemical  kinetic  effects. 
The  advantage  of  this  approach  compared  to  the  C/D  model  is  that  it  eliminates 
the  computer  generated  random  selection  procedure  and  is  therefore  more  efficient 
in  the  use  of  computer  time. 


Tracer  stimulus  response  techniques  have  been  previously  applied  to  the 

7  ,  8.9. 

study  of  chemical  reactor  design,  residence  time  distribution  in  combustors, 

mixing, ^  turbulence,^  etc.  However,  the  present  study  has  evolved  largely  from 

methods  developed  for  process  control  and  controller  design  in  which  stimulus 

12-14. 

response  is  used  to  identify  transfer  functions  and  fit'stochastic  models. 


In  spite  of  the  diversity  of  the  fields  and  applications,  the  techniques  and 
the  relevant  theories  are  basically  the  same  and  a  good  introduction  to  the 
subject  is  given  by  Godfrey. ^ 

In  the  experiments,  the  system  is  disturbed  by  a  suitable  input  signal,  and 
the  resulting  system  response  together  with  the  input  signal  are  analysed  by  a 
variety  of  techniques  to  gain  valuable  information  about  the  dynamics  of  the 
system. 

Examples  of  forcing  functions  are:  pulses,  steps,  sinusoids,  random  (noise), 
and  pseudo-random  binary  sequences  (PRBS).  Although  in  principle  the  same 
information  can  be  gained  irrespective  of  the  form  of  the  input  signal,  the  PRBS 
is  a  much  sure  useful  and  powerful  signal,  and  it  is  well  established  and  widely 
used  in  the  field  of  process  control.  On  the  other  hand  step  and  pulse  functions 
are  the  dominant  forms  of  input  signals  which  have  been  used  in  combustion  and 
reactor  design  studies,  with  little  or  no  practical  use  of  the  PRBS. 

As  will  be  shown  later,  the  advantages  of  the  present  approach  are  as 
follows : 

a)  the  PRBS  is  on  the  average  a  steady  state  function  and  the  combustor  can 
therefore  be  operated  at  practically  steady  state. 

b)  the  amplitude  of  pulse  or  step  signals  would  have  to  be  relatively  large 
compared  to  the  PRBS  disturbance  and  this  can  result  in  system  overloading 
and  unquantified  non-linearities. 

c)  for  a  given  signal/noise  ratio  the  experimental  time  is  relatively  modest. 

d)  many  disturbance  frequencies  are  generated  for  one  easily  generated  code. 

e)  special  methods  described  below  can  be  used  to  identify  the  non-li^ar 
elements  of  the  system. 

f)  quantitative  estimates  can  be  made  of  the  turbulent  mixing  'noise' 
contribution  to  the  response. 

The  generation  of  a  PRBS  may  be  easily  accoipplished  with  a  shift  register 
connected  as  shown  in  Fig.  4,  The  number  of  stages,  sequence  period  and  feedback 
connections  are  shown  in  Table  I.  The  sequence  may  be  started  with  any  non-zero 

V 
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value  in  one  or  more  of  the  registers,  and  it  will  be  noted  that  m  registers 
provide  a  sequence  period  of  2m-l.  A  second  class  of  PRBS  sequence  based  on 
quadratic  residue  codes  is  described  in  Ref.  16. 


THEORY  (LINEAR  SYSTEMS) 

As  illustrated  in  Fig. 5,  a  linear  system  with  input  signal  i(t)  having  an 
impulse  response  h(t)  produces  a  response  signal  y(t)  which  together  with  the 
noise  due  for  example  to  turbulence  n(t)  results  in  a  measurable  output  signal 
Z(t)  thus 

Z(t)  -  y (t)  +  n(t)  (2) 

The  process  response,  y(t),  is  given  by  a  weighted  sum  of  inputs  which  have 
occurred  in  the  past.  Thus  the  past  inputs  are  multiplied  by  the  weighting 
function  (or  impulse  response)  h(t)  as  given  by  the  convolution  integral: 

r°° 

y(t)  -  I  h(t. )  i  (t-t.)dt  (3) 

J  1  1 

o 

In  practice,  inputs  occurring  at  times  in  the  past  greater  than  the  settling 
time  Ts  have  no  effect  on  the  present  output  so 
TS 

y (t)  -  f  h(t£)  i  (t-t1)dt1  (4) 

Jo 

Thus  from  Equns  2  and  4. 


Z(t) 


h(t^)  i  (t-c^dt^  +  n(t) 


(5) 


The  cross-correlation  function  between  the  two  signals  i(t)  and  Z(t)  is  given  by 


R.  (T)  =»  lira 
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£  h  fL 

*  rr 


i (t)  Z(t+x)dt 
-T 


(6) 


and  similarly  the  auto-correlation  function  of  a  signal  is: 


R;i(x)  -  1 


n 


lira  i  f  i 
T-*°°  L 


i(t)  i  (t+x)dt 


(7) 


The  auto-correlation  thus  represents  the  similarity  of  a  signal  with  a  time 

shifted  version  of  itself.  For  example  if  a  wide  hand  noise  is  compared  with  a 

time  shifted  verion  of  itself,  only  a  small  time  shift  is  required  to  destroy  the 

\ 

In  the  limiting  case  of  white  noise,  the 


similarity,  which  does  not  recur. 


the  auto-correlogram  is  a  Dirac  delta  function  at  zero  phase  shift.  Pseudo-random 
binary  signals  have  levels  ±  1  (so  that  correlation  only  involves  addition  and 
subtraction)  and  may  only  switch  level  in  a  p re-determined  sequence  at  time 
intervals  At.  The  sequence  repeats  after  N  intervals  (period  T  =•  NAt)  yielding 
an  auto-corre lelogram  which  is  similar  to  a  truly  random  signal  as  shown  in  Fig. 6. 
It  is  this  property  which  makes  a  PRBS  a  useful  signal  to  study  the  system  responst 
since  the  relationship  between  the  input-output  cross  correlation  function  and  the 


input  auto-correlation  function  is 

T 

*■  r* 


»iZ<T> 


»<V  Ru<T-ti>dti  *  Ri„w 


where  R.  (t)  represents  the  cross-correlation  between  the  input  and  noise 
in 

signal  n(t).  For  a  long  correlation  time  (depending  on  signal/noise  ratio  as 
discussed  below)  and  assuming  that  the  noise  signal  is  unaffected  by  the  input 
signal,  we  obtain  R^(t)  =  0.  Now  assuming  R^  could  be  represented  by  a  Dirac 
delta  function,  Eqn.  8,  would  become  the  convolution  of  the  impulse  response  with 
the  Dirac  function  (for  times  less  than  T  =  NAt) .  Now  a  useful  property  of  the 
Dira  delta  function  is  that  its  convolution  with  the  function  yields  the 
function,  i.e. 


<  <5(t-t  )  f(t)  dt  =  f(t  )  if  a  <  tD<  b 

4  .  (9) 

=0  if  a  b  interval  does  not  contain  t  . 

o 

It  therefore  follows  that  the  cross-correlation  between  input  and  output  signals 


would  approach  in  the  impulse  resnonse 


R.  (T)  -  JC.h  (T) 
1Z  1 


where  depends  on  the  energy  of  the  input  pulse. 

However,  Fig. 5.  shows  that  for  a  PRBS  with  a  period  greater  than  the  settling 


time  (Ts) 


one  auto-correlation  triangle  of  height  (1  ♦  j^)  appears  in  the 


period  together  with  a  d  c  level  of  -  i.  Provided  that  At  is  small  with  respect 
to  the  time  constants  of  the  process,  i.e.  h(t)  is  approximately  constant  over 


the  width  of  the  auto-correlation  triangle,  the  mul tiplication  and  integration 
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of  Eq.  8.  produces 


5 

Riz(T)  *  (1  +  h  At  h  (T)  -  I  |h  (t]L)  dtl  + 


R.  (Z) 
in 


(11) 


* 


where  the  second  term  is  a  d.c.  level,  with  the  integral/  h(t^)dt^  being  the 

Jo 

steady-state  gain  of  the  process.  Over  a  long  integration  period  the  correlation 
» 

R-n  (Z)  "*'0.  Thus  the  system  dynamics  can  be  conveniently  studied  with  the  aid 
of  equation  (11)  . 

The  results  of  PRIJS  experiments  may  be  expressed  either  as  the  weighting 
sequence  (impulse  response)  for  discrete  points 


Nl 

Z^  »  E  w  i,  plus  bias  plus  noise 
t  j  9  1-s 


(12) 


or  in  terms  of  the  Z  transform  (which  is  a  closed  form  representation  of  the 
weighting  sequence). 


Z 

t 


-1  n 

w  r  b .  w 
1  1 

n 

l  ♦la.  w 
1  1 


i 

t 


♦ 


P 

1  E  d. 

n 


n 


t 


(13) 


where  w  is  the  shift  operator  defined  as  w  y^  =  y ^  k  represents  the  time 

delay  in  the  system  and  n  is  the  order  of  the  system  transfer  function.  The 

objective  of  the  Z  transform  model  analysis  procedure  is  to  evaluate  the 

coefficients  and  order  of  the  polynomials  which  give  the  best  fit  between  the 

model  predictions  and  the  experimental  data. 

Amongst  the  generally  available  identification  techniques  including 

correlation  methods,  instrumental  variables,  maximum  likelihood,  etc.  the 

1 8 

generalized  least  square  (GLS) .algorithm  h <s  proved  to  be  statistically 
sufficient  and  is  linear  in  its  formulation.  The  algorithm  can  provide  unbiased 
estimates,  a  desirable  feature,  as  data  required  for  the  formulation  of  simple 
dynamic  models  is  almost  always  corrupted  by  correlated  noise  and  the  need  to 
evaluate  unbiased  estimates  of  the  process  model  and  of  the  noise  becomes 
apparent.  Additionally  means  exist  to  apply  various  diagnostic  checks. 
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Equation  13  can  be  re-expres3ed  (upon  removal  of  any  system  integrations 


by  data  differencing)  in  the  following  convenient  form 

-k  AD 

A2t  =*  w  Blt  +  “C  nt  (U) 

where  A,  B,  C  and  D  represent  the  polynomials  in  w  * .  Thus  unless  ADC  *=  1 

the  process  parameter  estimates  using  simple  least  squares  will  be  biased. 
Hence  the  GLS  aims  to  transform  the  term  ADC  *n  to  an  uncorrelated  sequence 
via  the  following  steps: 

i)  The  process  parameters  (a^,b^)  are  initially  estimated  using  an  ordinary 
least  squares  algorithm. 

ii)  The  residuals,  viz 

et  ”  ^t  “  w  1  ®£t  (15) 

(where  *  denotes  estimated  values)  are  analysed  to  be  subsequently 
transformed  by  autoregression. 

Fet  -  nt 

...  .  .  ^ 

iii)  The  process  input  and  output  are  filtered  with  the  autoregression  F 

to  yield: 

F  *  .F  A 

-*■  FZt  and  it  =•  Fi^  (16) 


.  F  P 

iv)  i^  and  Z  are  used  to  obtain  a  new  least  squares  fit  to  commence  another 


iteration  cycle  proceeding  from  step  ii. 

The  final  model  i3  expressed  in  the  form 

A(FZt)  =  w"k  B  (Fit)  +  nt  (17) 

The  validity  of  the  estimated  process  and  noise  models  given  by  the  GLS 
algorithm  have  to  be  established  by  diagnostics  checks  for  the  statistical 
properties  of  the  models.  Currently  available  diagnostic  algorithms  utilize 
model  order  tests  including  determinant  ratio  test,  F-ratio  test,  loss  function 
analysis,  pole-zero  cancellation  and  tests  for  independence.  In  addition. 
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auto-correlations  of  the  residuals  and  cross-correlations  of  the  input  and  the 

residuals  can  be  employed  to  reveal  the  adequacy  of  the  fitted  model. 

Tentative  values  for  the  time  delay  and  model  order  can  be  obtained  via 

the  impulse  response  evaluated  directly  by  cross-correlation  of  the  input  and 

output  for  white  inputs,  Eqn.  10.  An  alternative  mean3  of  extracting  the  same 

14 

information  is  to  use  the  dete rminunt rat  10  test  .  The  test  is  especially 
valuable  in  that  it  limits  the  number  of  possible  model  orders  and  associated 
time  delays  prior  to  parameter  estimation  if  the  signal  to  noise  ratio  is 
reasonably  high. 

Loss  function  analysis  and  the  associated  F-ratio  test  are  employed  once 
the  process  and  noise  models  have  been  derived  by  the  GLS  algorithm.  They 
are  primarily  concerned  with  the  behaviour  of  the  error  function  in  the 
neighbourhood  of  taken  to  correspond  to  the  actual  system  order. 

Additionally,  residual  manipulation  in  terms  of  auto-correlation  functions  provide 
a  critical  assessment  of  the  structure  and  order  of  the  estimated  model.  The 
auto-correlation  function  of  the  residuals  is  employed  to  yield  an  indication  of 
the  whiteness  of  the  residuals.  The  cross-correlation  function,  on  the  other 
hand,  establishes  the  statistical  independence  of  the  residuals  and  the  process 
input . 

The  design  of  experin^nts  and  the  interpretation  of  the  results  require 

particular  consideration  which  may  be  illustrated  by  reference  to  the  more 

familiar  step  test  Fig. 7.  Indeed  it  is  often  convenient  to  carry  out  a  step 

test  to  obtain  criteria  (T,,T  ,T  and  noise  variance)  for  the  design  of  the 

dps 

PRBS  experiment.  The  criteria  to  be  fixed  are  the  digit  interval  At,  the  period 
NAt  and  the  signal  amplitude  ia,  determined  as  follows: 

1.  Clearly  the  PRBS  time  interval  At  must  be  short  compared  with  the  system 
time  delay  k  if  precision  in  the  measurement  of  k  is  to  be  attained. 

2.  The  time  interval  At  should  also  be  short  compared  to  the  estimated  time 

constant  T  of  the  system  (Fig. 7.)  although  there  is  no  advantage  in  excessivly 
P 

short  intervals  and  they  would  result  in  a  decrease  of  the  signal  to  noise 


ratio  of  the  output  measurements. 

3.  The  period  N  T  should  be  about  1.25  to  1.5  times  the  settling  time  T  .  If  the 

3 

lag  T  is  a  first  order  time  constant,  then  T  %  5T  .  Longer  periods  may 
P  s  p 

introduce  problems  due  to  signal  drift  although  this  can  be  reduced  by  trend 
analysis. 

4.  The  amplitude  Ot  should  be  the  largest  which  does  not  cause  undue  disturbance 

to  the  system.  This  follows  since  the  first  two  terms  on  the  right  hand  side 

2 

of  equation  11  increase  as  a  whereas  R[n(Z)  only  increases  linearly  with  a. 

5.  The  sampling  interval  of  the  output  measurement  should  preferably  be 
synchronized  to  the  start  of  PRBS  steps.  More  than  one  sample  per  At  interval 
will  result  in  more  data,  which  can  be  used  either  to  digitally  fitter  the 
input,  thus  improving  the  effective  signal/noise  ratio,  or  the  resulting 
weighting  sequence  can  be  smoothed  since  it  will  consist  of  more  points. 

The  static  gain  G  of  the  system  can  be  estimated  from  the  sum  of  the 
weighting  sequence  ordinates  or  by  the  relation: 

G  =*  ?  b./(l  +  ?  a.)  (18) 

1  1  l  1 

However,  the  signal/noise  ratio  is  better  if  a  simple  step  test  is  used  to 

19 

obtain  G  rather  than  a  PRBS 

It  is  shown  in  Fig. 8.  from  reference  20  that  simple  step  response  tracer 
results  are  insensitive  to  factors  such  as  recirculation  ratio  unless  very 
precise  measurements  are  obtained  in  the  vicinity  of  the  roll-over  point  in  the 
response  curve.  This  region  corresponds  to  the  peak  region  of  the  impulse 
response  curve,  and  it  will  be  shown  below  that  this  is  precisely  the  zone 
where  PRBS  methods  give  detailed  results. 

EXPERIMENT  AND  RESULTS 

Due  to  our  limited  understanding  of  the  reactor  network  and  turbulence/kinetic 

interactions  in  gas  turbine  combustors,  an  experimental  technique  which  can  be 

applied  in  such  complex  three-dimensional  flow#  is  plainly  required  to  study 

the  phenotrena.  As  a  first  step,  the  use  of  tracers  in  a  water  model  has  been 
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investigated.  One  important  difference  between  water  models  (or  air  models) 
and  hot  combustion  flows  is  that  the  former  are  isothermal.  If  geometrical 
similarity  is  maintained  and  Reynolds  numbers  are  chosen  to  ensure  that  the 
model  flow  is  turbulent,  then  as  a  result  of  the  well  known  jet  similarity 
characteristics,  the  flow  patterns  in  the  water  model  and  prototype  will  be 
closely  related.  Using  "partial  modelling"  (21)  methods,  in  which  the 
diameters  of  the  holes  in  the  combustor  wall  are  increased  in  proportion  to  the 
square  root  of  the  cold  jet  to  the  hot  chamber  gas  density  ratio,  the  relative 
velocities  and  flow  patterns  are  even  more  closely  modelled.  The  model  is  not 
abi_  to  take  full  account  of  the  effect  of  heat  release  on  the  turbulence 
structure,  nor  the  effect  of  the  turbulent  fluctuations  on  the  local  heat  release 
rate.  Thus  the  micromixing  'noise*  measurements  will  not  be  valid,  however, 
deductions  concerning  the  reactor  network  should  be  close  approximations  to  the 
situation  in  a  hot  combustor.  Some  extensions  of  the  techniques  to  direct 
measurements  in  hot  combustion  systems  will  be  discussed  later. 

In  the  experiments  a  full  scale  water  model  of  a  research  combustor  of 
dimensions  shown  in  Fig. 9.  was  constructed.  Salt  solution  could  be  injected 
at  different  points  through  a  hypodermic  needle,  although  in  the  results  reported 
here  the  salt  was  always  injected  through  one  of  the  primary  holes  since  we  are 
particularly  interested  in  the  primary  region  of  the  combustor.  The  salt  flow 
was  modulated  by  a  solenoid  valve  controlled  by  an  LSI-11  computer,  and  the 
salt  concentration  was  measured  by  a  conductivity  micro-probe  which  could  be 
positioned  accurately  anywhere  in  the  chamber.  A  detailed  description  of  the 
experimental  set  up  and  measuring  system  is  given  in  Ref.  22. 

The  investigation  employed  PRBS  and  step  signals  respectively  as  the 
system  modes  of  excitation  for  parameter  estimation  purposes.  The  particular 
PRBS  adopted  uaed  a  bit  interval  At  of  1*1  sec  with  a  15  bit  period  fed  to  the 
solenoid  valve.  A  number  of  measuring  stations,  representative  of  the  distinct 
mixing  zones  in  the  combustor  were  chosen  (as  shown  in  Pig. 9.)  for  identification 
and  parameter  estimation. 


A  typical  sequence  of  operational  steps  in  the  analysis  is  as  follows: 

A  determinant  ratio  test  was  used  to  limit  the  number  of  possible  model  orders 
for  the  range  of  time  delays  suggested  by  the  deconvoluted  impulse  response. 

GLS  parameter  estimation  was  then  applied  for  varying  model  orders  and  time 
delays.  Typically  noise  model  orders  of  10-15  were  adopted  in  an  iterative 
scheme  of  5-10  steps.  Diagnostic  checks  comprising  mainly  the  evaluation 
of  the  a  c  f  of  the  residuals,  the  c  c  f  between  the  input  sequence  and  the 
residuals,  and  pole-zero  cancellation  were  applied  for  varying  model  orders 
and  time  delays.  Parameter  estimation  and  subsequent  validation  were  carried 
out  in  a  strict  iterative  manner.  Finally  the  impulse  response  corresponding 
to  the  estimated  model  was  computed  and  this  was  compared  with  the  original 
weighting  sequence  obtained  from  the  cross-correlation  of  the  input/output  data. 

The  results  obtained  using  step  changes  in  salt  concentration  are  illustrated 
in  Fig.  10.  The  sharp  cut-off  at  Station  2  within  the  primary  jet  is  clearly 
demonstrated.  The  recirculation  in  the  primary  zone  which  behaves  a3  a  group 
of  stirred  reactors  as  defined  above  is  clearly  illustrated  by  the  response 
at  Stations  1,  7  and  8.  The  time  constant  of  the  exponential  decay  together 
with  the  known  primary  flow  rate  of  0.075  £/s  indicates  an  effective  volume  for 
this  group  of  reactors  of  0.2  l  which  is  consistant  with  the  anticipated  value. 
This  clearly  illustrates  that  turbulent  exchange  and  recirculation  involves 
virtually  the  whole  volume  of  the  primary  zone,  rather  than  one  sixth  which 
could  be  inferred  by  the  cyclic  symmetry.  The  series  response  of  the  secondary 
zone  i3  indicated  by  the  results  at  station  12.  Correcting  the  measured  time 
constant  for  the  effect  of  the  primary  zone  and  the  additional  secondary  flow 
gives  the  effective  volume  of  the  secondary  stirred  reactor  as  0.15  l. 

Turning  next  to  the  PRDS  technique,  some  of  the  Z  transform  results 
obtained  are  given  as  follows: 
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Station 

Location 

Order 

Time  Delay 

31 

a2 

bl 

b2 

7 

Primary 

1 

2 

-0.572 

- 

112.9 

- 

13 

Primary 

2 

l 

-0.753 

-0.117 

3.16 

-116.5 

9 

Secondary 

2 

3 

-0.293 

-0.293 

-78.3 

-61.8 

14 

Dilution 

2 

4 

-0.339 

-0.528 

-17.1 

-112.3 

The  corresponding  weighting  sequences  (impulse  responses)  are  shown  in  Fig.  1! 
Integrating  these  curves  gives  a  representation  of  the  step  response  which  can 
be  used  to  determine  the  effective  network  volumes  as  discussed  above.  However, 
the  integral  tends  to  mask  the  detail  which  can  be  seen  near  the  peak  of  the 
impulse  curves  as  clearly  illustrated  by  comparing  the  two  cases  at  Station 
14  (Fig.  11),  thus  it  appears  that  the  anticipated  advantages  of  the  PRBS 
method  may  be  realized  in  practice. 

The  experimentally  determined  impulse  response  at  Station  9  is  particularly 

interesting  since  this  point  is  on  the  opposite  side  of  the  axis  to  the  trace 

injection.  The  response  at  this  point  is  therefore  due  to  turbulent  exchange 

rather  than  convection,  and  the  fact  that  it  is  almost  identical  to  that  above 

the  axis  confirms  the  results  discussed  above  concerning  the  effective  volume 

of  the  primary  zone.  The  effect  of  this  recycling  reactor  group  size  is  to 

increase  the  throughput  at  blow-off,  that  is,  enlarge  the  stability  loop  and 

1 

altitude  performance,  as  discussed  in  Ref.  23. 

The  relationship  between  the  Z  transform  formulation  and  the  more  familiar 

l 

(at  least  tq  combustion  engineers)  time  constant  concept  is  as  follows  (24). 

1)  For  a  simple  first  order  system  with  no  delay,  the  continuous  system 


(ii)  For  a  delayed  first  order  system,  the  equation 
(1  ♦  Tp)y(t)  -  Gi(t-b-c) 


(21) 


corresponds  to 

(1  -dw  1)yt  -  (ojq  -  w1  w  l)it_b_1+  (22) 

where  =  G(l  -  6*  C)  and  =»  G(6  -  6^  C) 

(iii)  Similarly,  a  continuous  second  order  system  with  delay  obeying  the 
equation 

(1  -*•  T1D)(1  ♦  T2D)y(t)  -  i(t-b-c) 

where  and  T2  are  real  (rather  than  complex),  is  co-incident  with  the  discrete 
system  satisfying 

(1-S1w_1)(l-S2w"1)yt  -  (wo  -  c^w"1  -  V_2)it-b-U  (23> 

where  SL  -  e~1/Tl  ,  S2  -  e~1/T2  and 
0)o  -  G(T1  -  T2)“L  (T^l-Sj1-0)  -  T2(l-S2l_C)} 

Wj  -  G(Tl  -  T2)“1'  {(St  +  S2) (Tj  -  T2)  +  T2S21_C(l  +  Sx)  -  T1S11'C(1  +  S2)} 
w2  -  G  S1S2<T1  -  T2)_1'  {T2(l  -  S~C)  -  Tjd  -  S"C)} 

NB  b  is  the  number  of  whole  increments  At  and  c  is  the  fraction,  hence  b  ♦  c  »  del 
Thus  using  equation  (22)  the  Z  transform  coefficients  at  Station  7  are 
equivalent  to  the  continuous  function  Eqn.  (21)  with  <$  *  0.572,  c  »  0,  b  •  2 
gain  ■  -  263,  and  the  corresponding  time  constant  T  ■  (-  tn  6)  ^  *  1.79  secs. 

It  is  apparent  from  the  results  at  other  Stations  that  the  power  of  the  PRBS 
to  identify  two  or  more  time  constants  has  been  clearly  demonstrated. 

In  principle,  the  value  of  the  mixing  parameter  0(or  Tp)  can  be  crudely 
obtained  from  the  location  of  the  maximum  of  the  variance  of  the  step  response 
(see  Ref.  25),  however  in  well  stirred  systems  such  as  the  gas  turbine  combustor, 
the  location  of  the  maximum  is  so  close  to  zero  that  it  cannot  be  clearly 
distinguished.  The  equation  relating  the  location  of  this  maximum  to  3  for  a 
given  residence  time  Tg(-  l/a)  is  (25): 

t'  ,  ,  ^  -  (1/(0  -a)  )  £n  (  (0  +  a)/2d)  v 

i«y(t)  »max  % 


Fig.  12  shows  the  variance  of  a  sequence  of  concentration  steps  applied  to  the 
water  model  and  it  can  be  seen  that  the  maximum  variance  is  located  close  to  zero. 


This  is  consistent  with  a  high  value  of  6  as  follows:  Taking  a  typical 
residence  time  Tg  arfl  sec,  and  8  >  100  gives  the  time  of  the  maximum  variance 
1 1  <<  y (t) >>  <0.04  secs.  Fortunately  in  highly  stirred  systems  such  as  gas 

turbine  combustors  the  results  are  insensitive  to  the  exact  value  of  T^(or  8). 

In  fact,  the  combustion  efficiency  is  more  significantly  influenced  by  the 
uniformity  of  the  fuel  distribution,  whereas  the  production  of  N0^  may  be  more 
influenced  by  the  joint  probability  distribution  of  concentrations  and 
temperature  (which  is  at  present  unknown),  than  by  the  exact  value  of  t^. 

Since  probability  distributions  tend  towards  gaussian  when  they  depend  on  the 
combination  of  several  other  probability  distributions,  almost  whatever  the 
shape  of  the  separate  distributions,  it  is  likely  that  the  effects  of  the  pdf 
on  pollutant  formation  can  be  represented  by  a  set  of  about  five  parallel  stirred 
reactor  loops,  the  size  and  operating  conditions  of  each  weighted  according  to 
a  portion  of  the  gaussian  curve. 

If  an  accurate  estimate  of  Tp  throughout  the  combustor  is  required  it 

should  be  possible  to  extract  this  from  the  noise  term  in  Eqn.  (13).  Although 

the  auto-regressive  model  F  used  in  the  CLS  algorithm,  Eqn.  (16),  is  often  of 

“I  -1 

high  order,  this  is  to  ensure  that  F  *  (ADC  )  provides  an  adequate 
representation  of  the  rational  transfer  function  (ADC  .  Estimates  of  the 
noise  model  coefficients  can  however  be  obtained  by  including  an 

additional  step  in  the  GLS  algorithm.  This  would  consist  of  computing  the 
least  squares  estimate  of  d^,  in  the  model 


e 


t 


D 

C  nt 


(22) 


where  the  residuals  nt  are  considered  as  input  variables  and  the  deterministic 

-k  -l 

production  errors  e  ■  zfc  -  w  A  Bi^  as  output  variables.  Both  and  nt 
are  generated  within  the  GLS  algorithm  and  hence  could  be  derived  from  the 
estimates  associated  with  Eqn.  (22). 

v 
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Further  developments  of  tracer  techniques  currently  being  investigated 
at  Sheffield  are  based  on  hot  combustion  systems.  In  the  first,  a  rapid 
response  (x»*0.1  sec)  infra  red  CO^  analyser  has  been  developed  to  study  the 
overall  chemical  reaction  dynamics  in  a  coal  fired  fluidized  bed  combustor. 

The  CO^  analyser  is  used  to  monitor  the  conversion  achieved  in  the  bed  in 
response  to  a  PRBS  input  of  coal  or  various  reacting  gases.  In  the  second 
study,  the  flow  reactor  network  in  a  hot,  high  intensity  combustor  is  to  be 
investigated  using  a  mercury  vapour  tracer  as  described  by  Topps  (26).  He 
demonstrated  that  the  vapour  may  be  generated  by  passing  brief  electrical 
discharges  through  twin  bore  refractory  tubing  filled  with  mercury  amalgam, 
and  subsequently  detected  by  the  absorbtion  of  light  generated  by  a  hollow 
cathode  mercury-neon  lanp.  In  both  studies,  since  only  the  a  c  response 
signals  are  required,  long  term  stability  is  relatively  unimportant.  The 
experimental  results  using  the  methods  detailed  in  this  paper  will  be  reported  in  d 
course. 

At  this  point,  the  application  of  the  PRBS  stimulus  technique  to  the 
measurement  of  the  acoustic  admittance  (response  function)  of  burners,  rocket 
motors  and  propellants,  boiler  and  furnace  chambers,  etc.  should  be  mentioned. 

The  particular  advantage  is,  of  course,  associated  with  the  large  number  of 
disturbance  frequencies  generated  by  the  single  code. 

THEORY  (NON-LINEAR  SYSTEMS) 

The  representation  of  a  mixing  process  followed  by  the  non-linear  Arrhenius 

rate  equation  in  cascade  with  the  instrumentation  network  is  difficult  to 

analyse  experimentally  because  of  the  complexity  of  the  system.  However,  recent 

27  28 

results  in  non-linear  identification  theory  ’  have  shown  that  systems  of 
this  structure  can  be  identified  in  terms  of  the  individual  component  subsystems 
from  measurements  of  the  input  and  noise  corrupted  output  signals  only. 

Classically  identif ication  techniques  for  non-linear  systems  are  usually 
based  on  the  Volterra  series  description: 
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Ct)  g1(T1)  i(t-T1)dT1  +JJ  g2(T1,T2)  i(t-T1)  i  (t-T2)dT1  dT, 


+ . +/  .  .  . /gm  (Tx,  .  .  .  Tn)  i  (t-T1)  .  .  . 


•  •  •  •  i  (t  -Tm)dTl  •  •  dTm  +  n(t) 


E  w,  (-t)  +  n (t) 
k=l 


(23) 


.  20  .  _ 
or  the  related  Wiener  expansion  ,  where  the  function  g:  T,  .T,.  .  T  )  is  termed 

1  /  n 

the  Volterra  kemal  of  order  m  .  The  first  term  in  the  Volterra  series  is  the 
well  known  convolution  integral  associated  with  linear  systems  and  the  higher 
order  terms  reflect  the  non-linearity  of  the  process.  Unfortunately,  identification 
of  the  Volterra  kernals  often  involves  an  excessive  computational  requirerrent  and 
the  physical  interpretation  of  the  final  model  in  terms  of  the  components  of  the 
original  system  is  usually  difficult  to  achieve.  Consequently,  very  few 

authors  have  considered  the  identification  of  practical  systems  using  these 
techniques. 

27  28 

However,  recent  results  ’  have  shown  that  systems  composed  of  linear 
dynamic  and  static  non-linear  elements  can  be  identified  using  relatively  simple 
extensions  of  well  established  linear  techniques  in  a  manner  which  preserves 
the  system  structure  and  provides  estimates  of  the  individual  component  subsystems. 

In  the  present  study  only  the  cascade  non-linear  system~Tl  lustrated,  in  Fig. 13. 
which  has  the  structure  of  a  stirred  reactor  network  will  be  considered.  The 
identification  problem  can  be  defined  as  identification  of  hj(t)  the  mixing 
dynamics,  F[*j  the  Arrhenius  rate  equation  and  h2(t)  the  instrumentation  network 
from  measurements  of  i^(t)  and  z^(t). 

If  the  input  i*(t)  is  a  separable  process,  and  separability  is  preserved 
under  linear  transformation  for  i^(t),  then  the  first  degree  cross-correlation 
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function  is  given  by  * 

R,  .  (e)  =  CFG(f  h2(0)h1(t1)  (£-9  -  T1>  d9  dT1  +  ^E)  (24) 

Similarly,  correlation  between  the  input  squared  and  the  output  yields  the 

second  degree  correlation  function,  when  separability  is  preserved  under 

27  28 . 

double  non-linear  transformation.  ’ 


R-2?  (e)  =  CFFG[|[  h2(9>h1(T1)  hl  (x2)  i1(t-0-T1)  ix  (t-0-  x2) 


i?  (t-£)  dT.dT„d0  +  R.2  (e) 
1  12  i»n 


(25) 


For  the  case  when  i^(t)  =  i(t)  +  b  where  i(t)  is  a  zero  mean  gaussian  white 
process  and  b  is  a  non  zero  mean  .evel  equations  (24)  and  (25)  reduce  to 

^  zj  (e)  =  CFGJ  h1(T1)h2(e-T,)dx1  (26) 

Ri2z j  <e>  “  CFFcfh2(Tl)hl  (e~VdTl  (27> 

CFG  -  +  2Y2bJh1(0)d0  +  3Y3fhf  (0)  +  3Y3 b2  f  h1(T1)h1(x2)dx1dT2  +  .  .  (28) 

CFFG  =  2Y2  +  by2bf  hl(9)  dQ  +  •  •  •  (29) 

where,  provided  h  ( t )  is  stable  bounded-inputs-  bounded-outputs  C  and  C„_„  are 
1  ^  tO  FFG 

constants,  F[*]  =  E  y.(-)1  >  R.  (e),  R?  n(£)  tend  to  zero  when  n(t)  is 

i=l  ln  1 

independent  of  the  input  and  the  superscript  '  is  used  to  indicate  a  zero  mean 
process . 

/  Of.)  ( 0  7) 

The  estimates  of  equations''  and  are  quite  independent  of  the  non¬ 

linear  element  F[*]»  the  Arrhenius  rate  equation  in  this  case,  except  for  the 
constant  scale  factors  C  and  C  Correlation  analysis  thus  effective!'’ 

decouples  the  identification  problem  into  two  distinct  steps;  identification 
of  the  linear  subsystems  and  characterisation  of  the  non-linear  element. 

Since  we  effectively  have  two  equations  in  two  unknowns,  estimates  of  the 
individual  linear  subsystems  y^h^it),  ^^(t),  where  y^  and  y2  are  constants* 
can  be  obtained  from  equations  (26)  and  (27)  using  a  least  squares  decomposition 
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ilgorithn??  Once  the  linear  subsystems  have  been  identified  the 


input  and 


output  of  the  non-linear  element  are  effectively  known.  The  problem  is  therefore 

reduced  to  fitting  an  appropriate  function  to  the  non-linear  element  by 

minimising  the  sum  of  squared  errors. 

The  results  of  eqn.  (26)  and  (27)  also  provide  information  regarding  the 

system  structure,  or  position  of  the  non-linear  element  with  respect  to  the 

linear  subsystems.  If  the  system  is  linear,  =  0,  k  ^  1  it  can  readily  be 

shown  that  R.2  ,  (e)  =  0  H  e  and  this  provides  a  very  convenient  measure  of 
l  z1 

linearity.  If  the  1st.  and  2nd.  degree  correlation  functions  are  equal  except 
for  a  constant  of  proportionality,  then  the  system  must  have  the  structure  of 
the  Hammerstein  model  with  h^(t)  =  <S(t).  However,  if  the  2nd.  degree  correlation 
function  is  the  square  of  the  1st.  degree  correlation  function,  except  for  a 
constant  scale  factor,  the  system  must  have  the  structure  of  the  Wiener  model 
with  1^(0  =  6(t).  Finally  if  none  of  the  above  conditions  hold  the  system 
may  have  the  structure  of  the  general  model  but  this  must  be  confirmed  by  residual 

30 

analysis  or  application  of  an  algorithm  by  Douce  . 

Although  the  results  presented  above  are  relatively  simple  to  compute, 
limitations  of  the  input  transducers  may  dictate  that  the  input  should  be  a 
pseudorandom  sequence  rather  than  a  gaussian  white  noise  process.  Unfortunately, 
the  results  of  eqns  (26)  and  (27)  do  not  hold  for  pseudorandom  input  sequences 
because  of  anomalies'®^ associated  with  the  higher  order  correlation  functions  of 
these  inputs.  However,  by  injecting  a  compound  input  defined  as  i^(t)  = 

V 

X^(t)  +  *2(0  analogous  results  can  be  obtained  by  isolating  the  correlation 

functions  associated  with  the  outputs  Wj,(t)  and  W2(t)  eqn.  (23),  of  the  fi:^t 

31 

two  Volterra  kernels  to  yield 


»;  <0  -  BjYjhjU-s)  h2(0)  de 
Rx;x'  w'  (£)  '  26l62Y2.f  hl  <s"0)  h2<e>de 


where  X ^  ( t )  and  XjU)  are  uncorrelated  Rx  x  ^  =  0  H  X,  zero  mean  pseudorandom 


19  - 


sequences  with  autocorrelation  functions  x  (X)  «  f^6(X),i  -  1,2.  The 

i  1 

disadvantage  of  this  approach  is  the  requirement  to  use  multilevel  compound 

inputs  K.  i,  (t)  in  order  to  isolate  IL-i  ,  (£)  and  ILW,  ,  (e)  . 
i  1  *lwl  X1*2W2 

Since  the  results  of  equations  (30)  and  (31)  are  dependent  upon  X^(t)  and  X2(t) 
having  a  zero  mean  value  an  obvious  choice  of  input  would  be  a  compound  ternary 
sequence.  Pseudorandom  binary  sequences  can  however  be  used  but  the  non-zero 
mean  level  of  these  sequences  introduces  a  small  bias  into  the  estimate  of  eqn .  (31) 
This  bias  tends  to  zero  as  the  sequence  lengths  are  increased  and  will  be 
negligible  in  most  applications. 

Thus  by  applying  the  above  results  to  the  stirred  reactor  network  it  should 
be  possible  to  identify  the  mixing  dynamics,  parameters  associated  with  the 
Arrhenius  rate  equation  and  the  dynamics  of  the  instrumentation  network  and  this 
would  undoubtedly  provide  a  considerable  insight  into  the  operation  of  the  process. 

CONCLUSION 

The  problem  of  macro- and  micro-mixing  in  combustors  has  been  studied  by 
PRPS  stimulus  tracer  response  techniques  using  a  water  model.  The  potential 
advantages  of  the  method  have  been  clearly  demonstrated  and  future  applications 
to  the  measurement  of  non-linear  system  parameters  presented. 
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